%% pde_MF_s
clc;
clear;

index = 1;

syms x y nu nu_r C0 Ca Cd
switch index
    case 1
        u1 = -sin(pi*x).^2.*sin(pi*y).*cos(pi*y);
        u2 = sin(pi*x).*cos(pi*x).*sin(pi*y).^2;
        w3 = u1+u2;
        p = sin(pi*x).*cos(pi*y);
    case 2
        u1 = pi*sin(pi*x).^2.*sin(2*pi*y);
        u2 = pi*sin(2*pi*x).*sin(pi*y).^2;
        w3 = u1-u2;
        p = 10*cos(pi*x).*sin(pi*y);
    case 3
        u1 = 10*x.^2.*(x-1).^2.*y.*(y-1).*(2*y-1);
        u2 = 10*y.^2.*(y-1).^2.*x.*(1-x).*(2*x-1);
        w3 = u1-u2;
        p = 10*(2*x-1).*(2*y-1);
end

u1dx = diff(u1,x);
u1dy = diff(u1,y);
u2dx = diff(u2,x);
u2dy = diff(u2,y);
w3dx = diff(w3,x);
w3dy = diff(w3,y);
pdx = diff(p,x);
pdy = diff(p,y);
u1dxdx = diff(u1dx,x);
u1dydy = diff(u1dy,y);
u2dxdx = diff(u2dx,x);
u2dydy = diff(u2dy,y);
w3dxdx = diff(w3dx,x);
w3dydy = diff(w3dy,y);

f11 = -(nu+nu_r)*(u1dxdx+u1dydy) + (u1*u1dx+u2*u1dy) + pdx - 2*nu_r*w3dy;
f12 = -(nu+nu_r)*(u2dxdx+u2dydy) + (u1*u2dx+u2*u2dy) + pdy - 2*nu_r*(-w3dx);
f2 = -(Ca+Cd)*(w3dxdx+w3dydy) + (u1*w3dx+u2*w3dy) + 4*nu_r*w3 - 2*nu_r*(u2dx-u1dy);

disp("===u1导数===");
matlabFunction(simplify(u1dx), "Vars", {x,y})
matlabFunction(simplify(u1dy), "Vars", {x,y})
disp("===u2导数===");
matlabFunction(simplify(u2dx), "Vars", {x,y})
matlabFunction(simplify(u2dy), "Vars", {x,y})
disp("===w3导数===");
matlabFunction(simplify(w3dx), "Vars", {x,y})
matlabFunction(simplify(w3dy), "Vars", {x,y})
disp("===p导数===");
matlabFunction(simplify(pdx), "Vars", {x,y})
matlabFunction(simplify(pdy), "Vars", {x,y})
disp("===定常MF方程===");
matlabFunction(simplify(f11), "Vars", {x y nu nu_r C0 Ca Cd})
matlabFunction(simplify(f12), "Vars", {x y nu nu_r C0 Ca Cd})
matlabFunction(simplify(f2), "Vars", {x y nu nu_r C0 Ca Cd})
